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Abstract 

We introduce an adaptive regularization approach. In contrast to conventional 
Tikhonov regularization, which specifies a fixed regularization operator, we esti- 
mate it simultaneously with parameters. From a Bayesian perspective we estimate 
the prior distribution on parameters assuming that it is close to some given model 
distribution. We constrain the prior distribution to be a Gauss-Markov random 
field (GMRF), which allows us to solve for the prior distribution analytically and 
provides a fast optimization algorithm. We apply our approach to non-rigid image 
registration to estimate the spatial transformation between two images. Our eval- 
uation shows that the adaptive regularization approach significantly outperforms 
standard variational methods. 



1 Introduction 

Tikhonov regularization has been a standard tool to tackle ill -posed problems |[Tl|2l. One often 
minimizes an objective function regularized with a smoothness constraint. 

£;(u) = i:i(0,u) +w||Pu||^ (1) 

where D{0,u.) is a measure of how well the solution u fits the given data O and ||Pu||^ is a 
regularization term that penalizes some properties of u (e.g. lack of smoothness, when P is a 
derivative operator). Parameter w is a trade-off between data fitness and regularization. General 
non-quadratic forms of the regularization term have also been used. 

Optimization of the regularized objective function is often challenging, because the algorithms tend 
to get stuck in local minima. To overcome this problem one can adjust the value of the regularization 
parameter w. Large values of w allow to overcome some local minima, but result in overconstrained 
solutions. Thus, w has to be chosen big enough to avoid local minima, but small enough to allow 
flexibility on u. Multiple strategies to select w have been proposed, including various heuristics, 
slow annealing, cross validation and Bayesian estimation |(2ji3J. In many cases, there may not be a 
single w adequate to achieve a reasonable solution. 

Instead of searching for an optimal w with a fixed regularization operator, we estimate the reg- 
ularization operator P, treating it as an unknown parameter We first consider the regularization 
framework from a Bayesian perspective, where the regularization term comes from a Gaussian prior 
on u, and P^P is the inverse covariance matrix (or potential matrix). Instead of fixing the prior 
distribution (equivalent to fixing P), we estimate it assuming that it is close to some given model 
distribution. As we shall show, this allows flexibility on the prior distribution and leads to adaptive 
regularization. We constrain the prior distribution to be a Gauss-Markov random field (GMRF) due 
to the following reasons: a) a GMRF on a finite lattice has a shift invariant covariance, which allows 
us to solve for the covariance matrix analytically, b) the shift invariance property is consistent with 
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derivative based regularization, and c) the known eigenstmcture of the covariance matrix allows fast 
optimization. 

We introduce our new regularization approach from a general Bayesian perspective and then con- 
sider in detail the specific problem of non-rigid image registration. In non-rigid image registration 
one needs to find a non-rigid transformation that aligns two given images. Non-rigid image regis- 
tration is one of the key problems in computer vision with multiple application including motion 
correction, cross modality image fusion and atlas constraction |4, 5|. The rest of the paper is orga- 
nized as follows. In Sec|2]we define a general adaptive regularization framework and describe the 
fast optimization algorithm. In Sec. [3] we overview the non-rigid image registration problem and 
show how to apply the adaptive regularization approach. In Sec.|4]we evaluate our algorithm. In 
Sec.|5]we compare our algorithm to related methods. Sec.|6]concludes the work with discussions. 



2 Method 



2.1 Bayesian formulation 

From a Bayesian perspective the regularization approach is equivalent to the maximum a posteriori 
(MAP) estimation, that is to maximize 

maxp(u|0) cx p(0|u)p(u) (2) 

or equivalently to minimize the following objective function 

min£'(u) — — logp(0|u) — logp(u) (3) 

where the first term (the negative log-likelihood) is the error function D{0, u) and p(u) is a prior 
distribution on u. In case of the quadratic regularization term (our case), the prior p(u) is a Gaussian 
distribution 

p(u) = ^ =e-^"^^"'" (4) 

v/(27r)^det(I]) 

Defining the inverse covariance matrix (also called potential matrix) as = P^P and substitut- 
ing p(u) in Eq.[3] we achieve 

min£;(u,P) = -D{0,n) + \ |lPu||' - 1 logdet(P^P) + log(2^) (5) 
W z z z 

The last two terms are constants if P is fixed. The weight w here comes from the likelihood function 
(this is equivalent to the regularization weight in Eq.[T]). We prefer to use the covariance form instead 
of the operator form, and rewrite Eq.|5]as 

min£:(u,i;"i) = -DfO. u) + ^u^E^^u - i log detfS^^) + const (6) 
w ' 2 2 

So far we have just reformulated the regularization approach from a Bayesian perspective without 
any modifications. Now, instead of assuming a fixed inverse covariance matrix Yr^ (equivalent to 
assuming a fixed operator P), we shall estimate it. 



2.2 Adaptive regularization 

We assume that the prior distribution p(u) is unknown, but is close to the model distribution q{\\) in 
terms of Kullback-Leibler (KL) divergence. We rewrite the MAP problem as 

mini?(u,p) = -logp(0|u) - logp(u) ^ KL(j)\\q) (1) 

where KL{p\\q) is the KL-divergence between the unknown prior distribution p and the model 
distribution q 

KL{p\\q) = I p(u)log^du (8) 

We want to minimize Eq. [T] simultaneously with respect to u and p. The prior distribution p is 
multidimensional in general. To simplify and stay consistent with conventional regularization, we 
assume that p and q are zero-mean multivariate Gaussian with covariances E and il respectively: 
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p ^ A/'(0, E), g ~ A/'(0, 51). One advantage of such assumption is that the KL-divergence has an 
analytical form 

KL{p\\q) = i (tr(f7-iE) - logdet(O-i) + logdet(E-i) - N) (9) 
which we substitute in Eq.|7]to obtain 

min£;(u, = -D{0, u) + iu^E^^u + i + const (10) 

w 2 2 

We only need to find E to fully specify the prior distribution p. We can analytically solve for E 

in at least two cases: a) when both E and ft are diagonal, b) when E and ft commute (share the 

eigenvectors). Indeed the second case arises naturally in regularization theory. 



2.3 Markov random field and shift invariance 



We constrain the prior distribution to be a Gauss-Markov random field (GMRF), which implies 
that the inverse covariance matrix is shift invariant and is diagonalized by the Fourier basis on the 
discrete lattice |6|. Matrices E^^ and 51^^ are symmetric positive definite by definition and allow 
spectral decomposition 

E^i^QAQ^, n^^QKQ^ (11) 

where A = d[Ai, .., Ajv], VA„ > is a diagonal matrix of the eigenvalues of E^^ and K = 
d[fci, .., /cjv], Vfc„ > 0, is a diagonal matrix of the eigenvalues of 51^^. For the shift invariant matrix, 
the eigenvectors Q = [qi, ••, qjv] have a known form. Depending on the boundary conditions, Q 
is a discrete Fourier basis (circular boundary conditions) or discrete cosine transform (DCT) basis 
(Neumann boundary conditions) (l6]|7l. The assumption of the GMRF structure has multiple advan- 
tages. First, estimation of the inverse covariance matrix simplifies to estimation of its eigenvalues. 
Second, the product Q^x is a multidimensional DFT (or DCT) transform and can be computed 
fast in 0{N log N). Finally, the shift invariant structure of the covariance matrix is consistent with 
standard derivative-based regularization operators. 

Common regularization operators P, such as first, second, and higher order derivative operators are 
shift invariant by definition. In the discrete case, this means that the matrix 

E"i=P'^P (12) 

is a Toeplitz plus nearly Hankel matrix fV,^. Such a highly structured matrix is known to be diago- 
nalized by a Fourier basis, and the estimation of P^P reduces to an estimation of its eigenvalues L9J . 
Substituting Eq.fTTIinto Eq.[TOl we obtain 

min£;(u, A) = ^0(0, u) + ^ ^ A,(qf u)^ + ^ E | (1^) 
We equate the gradient of the function to zero, and solve for A to obtain; 

A. = ^ (14) 

|q/ u| 

where | • | denotes the absolute value. The solution for the eigenvalues A is guaranteed to be positive. 
Substituting A back into Eq. [T3] we obtain 

minS(u) = -15(0. u) +VyfcI|qfu| (15) 
w ^-^ 

At this point, we have analytically solved for A (which uniquely specifies the prior distribution p) 
and eliminated it from the objective function. The final form of the objective function includes the 
regularizer that penalizes the absolute value of Fourier (or DCT) coefficients of u weighted by K. 
We still need to specify the eigenvalues K of the model distribution q. 

To choose the model distribution q, and in particular its inverse covariance matrix 51^^, we follow 
a standard regularization approach. For instance, Laplacian ||Au||^ is a popular regularizer, which 
suggests to use 17 ~^ = as the inverse covariance matrix of the model distribution. The eigenval- 
ues of the discrete Laplacian (squared) on a regular grid (with Neumann boundary condition |i7J) in 

'„1 iT 



ID are = \k\, ..kj^f, fc„ = 2(1 - cos{TT{n - l)/N)), n = l,..,N. 



By shift invariant matrix, we denote a Toeplitz plus nearly Hankel matrix (due to the boundary conditions). 
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2.4 Optimization 



To optimize the objective function in Eq.[T5] we take advantage of the fact that the original objective 
function (Eq.fToli had a quadratic regularization term, and split the optimization into two steps: 

• findE^^: A = K^/^ ^m.g(|Q'^u|)"\ = QAQ'^, (16) 

• minimize £;(u) = D{0, u) + wu^E^^u (17) 

The first step, solution for E^^, has a closed form, whereas the second step requires iterative mini- 
mization (unless D{0, u) is quadratic). We shall briefly state one of the standard fast optimization 
approaches to minimize the function with quadratic penalty term. 

We equate the gradient of the objective function in Eq.[l7]to zero 

VE{u) = VD{0,u) +wY.-^u = (18) 

Note that the gradient consists of the non-linear part {WD{0, u)) and the linear part (E^^u). We 
artificially introduce a time-step derivative to the right-hand side of Eq.fTSjas 

VL»(0, u*) + wE~iu*+i = -(u*+i - u*)/7 (19) 

which converges to zero in equilibrium. Here 7 is a time step parameter (similar to the gradient 
descent step size parameter). Note that the linear part is at time t + 1, whereas the non-linear part is 
kept at t. Solving for u*+^, we achieve 

u*+i = Q(I + 7u.A)-iQ^(u' - jS/D{0, u*)) (20) 

Iterating on u we achieve faster minimization compared to the first-order minimization methods, 
with no need for Hessian computations or approximations. The eigenvector matrix Q is never con- 
structed explicitly. The matrix vector products Q^x and Qx are forward and inverse multidimen- 
sional DCTs respectively. Finally, combining Eq.[T6land Eq.|20]into a single step, we achieve our 
fast optimization algorithm by iterating on 

""■-Q miq'IT"'^) ^'''"'--'^"'"''-''' 

where d | • | denotes a diagonal matrix formed from the right-hand side vector. 



3 Non-rigid Image Registration 

Image registration is one of the key problems in computer vision. The goal of image registration is 
to find a spatial transformation that aligns two images. The most challenging cases of image reg- 
istration occur when the underlying transformation is non-rigid. Non-rigid image registration has a 
wide range of applications, including motion correction, change detection over time, cross modality 
image fusion, inter-subject comparison, atlas constraction, registration-based segmentation, motion 
estimation and tracking L4j|5l. 

As far as non-rigid transformation is a broad class of nonlinear transformations, non-rigid image 
registration is an ill-posed problem. To tackle the problem one can use either parametric or non- 
parametric (also called variational) approaches. Parametric image registration specifies a parametric 
model of the transformation (e.g., locally affine or B-splines), which explicitly constrains the trans- 
formation |10|. This, however, significantly limits the range of admissible transformations, and is 
not adequate for complex non-rigid transformations. Non-parametric image registration estimates 
the transforamtion as an unknown function using variational calculus. One often uses a regularizer 
on the transformation to make the problem well-posed. 

One of the standard non-parametric approach is to minimize the following objective function 01 111121 

E{u) = D{I,J^;u)+w\\Auf (22) 

where D{I, J^; u) is a similarity measure between the images / and J, and u is a displacement field 
that aligns J onto /. Regularization of the displacement function is often called the competitive reg- 
ularization, in contrast to the incremental regularization, which penalize increments (or evolution) 
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of the displacement function fT3|. Unfortunately, regularization of the displacement function (and 
its increments) do not guarantee a diffeomorphic transformation. A diffeomorphic transformation 
ensures that an inverse transformation and its derivatives exist and are smooth functions, which can 
be required in medical images. In a few cases when the standard regularization does not produce dif- 
feomorphic transformation, it can be post-smoothed to ensure the invertibility 1 14|. A more elegant 
approach to ensure diffeomorphism is to consider the transformation as a solution of the ordinary 
differential equation 1 15 , 16 1. Such diffeomorphic image registration methods explicitly account for 
invertible transformation, but suffer from large computational complexity. The detailed overview of 
the non-rigid transformation models in image registration is beyond the scope of this paper Here, 
we apply our regularization approach to estimate the displacement function, which is most similar 
in formulation to the competitive regularization in Eq.l22l 

Adaptive regularization: Following our adaptive regularization method, we minimize the follow- 
ing objective function 

mmE{u) = -D{I,J^;u)+y ^,\q[u\ (23) 

where !?(/, Ju;u) is a similarity measure between the images / and J. We shall 
use the sum-of-squared-differences (SSD) similarity measure, defined as £)(/, Ju;u) = 
(/(xi) — J(xj + u(xj)))^. The choice of the similarity measure is itself a research topic in 
image registration. Here, we use a simple similarity measure to concentrate primarily on the trans- 
formation estimation. Any other similarity measures can be easily applied for -D(/, Ju; u). We 
remind that are the DCT bases, which corresponds to the Neumann boundary condition on u. 
Such condition is often the most appropriate boundary condition for natural images. 

Implementation: To optimize the objective function in Eq. |23j we follow our iterative update 
according to Eq.lJT] To clarify the algorithm in multidimensional cases (e.g., for 3D images), we 
rewrite Eq.|2T|as 

|Q^u*| = ^DCT^ini) + DCT^{vlI) + DCT^{ui) + e 

u*+i = IDCT ( (|qtL^|'^"'{^k) ■ - ' * ^ ^' ' 

where DCT and IDCT denote forward and inverse multidimensional (in this case 3D) discrete 
cosine transforms, and all operations are elementwise. We add a small positive constant e (e.g. 
machine precision) to avoid devision by zero. Ua; , Uj, and are the 3D arrays of cooresponding 
voxel displacements. Array K denotes the eigenvalues of the multidimensional Laplacian with el- 
ements: hjp = IC{i,N^) + IC{j,Ny) + K:{p,N^), where /C(n, A^) = 2(1 - cos(^^^^)), and 
Nx, Ny, Nz are the image dimensions. Finally for the SSD similarity measure, the gradient is 
ViD{0, u*) = ( J(x + u) - /(x))Vi J(x + u), for i^x,y, z. 

4 Results 

We have implemented our algorithm in Matlab, and tested it on a AMD Opteron CPU 2GHz Linux 
machine with 4GB RAM. We used the BrainWeb Tl-weighted MRI images (180 x 216 x 180 
voxels) to test the algorithm [17J. We normalized image intensities to the [0,1] interval before 
registration. The stopping criterion was either when the algorithm reaches 1000 iterations or the 
objective function tolerance drops below 10^*. 

To simulate the synthetic spatial deformation we put a uniform grid of control points (with 15% 
spacing) over the image, randomly perturb them {JV{Q, a = 6)) and interpolate the image using thin 
plate spline [18 1. This way we obtain known smooth locally varying synthetic deformations. 

We compare our algorithm to the curvature based variational registration fTT, TSl, which has a 
quadratic Laplacian regularization term as in Eq. |22l To evaluate the registration performance, we 
compute the root mean squared error (RMSE) between the true and estimated transformations as 
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(a) reference (b) true transform (c) source (d) result 




(e) result (composite) (f) found transform (g) result (Laplacian) (h) transform (Laplacian) 



Figure 1 : Non-rigid registration example with complex synthetic transformation. We register the 
source image (c) onto the reference image (a). The registration result (d) is accurate. Figure (e) 
shows the composite view of registered image and the contours extracted form the reference image. 
The estimated transformation (f) is very similar to the true transformation (b) with Ermse — 0.52 
transformation error. Note, that we do not include the area outside the scull into the error com- 
putation to avoid the boundary error influence. We compare our results to the Laplacian-based 
regularization (g,h), which has less satisfactory performance with Ermse = 2.21. 



£ RAISE ^ jnJ2 \\^true - Uestimated\\'^ ■ We do not includc the area outside the scull (roughly 
found by thresholding) for the RAISE computation to avoid the boundary error influence. 

First, we demonstrate the registration performance on 2D images (216 x 180). Figure[na,b,c) shows 
the reference image, its synthetically deformed version (source image) and the true transformation 
(which gives 7.32 initial RMSE). Figure [nd,e,f) shows our registration result. Our algorithm ac- 
curately estimates the transformation with ermse = 0.52. Laplacian-based registration (g,h) pro- 
duces less satisfactory result with Eumse — 2.21. Notice, that the Laplacian-based approach did 
aUgn some less challenging parts of the image, but got stuck in a local minima. 

Figure |2] demonstrates the full volume 3D non-rigid image registration. The registration perfor- 
mance is accurate. We conducted several tests to evaluate and compare our approach to the standard 
quadratic (Laplacian-based) regularization. We performed 100 non-rigid 3D image registration ex- 
periments with random non-rigid deformation (similar to the ones in Fig. [TJ)) initialization at each 
run. The average initial transformation RMSE was 8.36. We did 20 experiments for each of the 
five different regularization weight values. We used w — [0.5, 1.2, 2, 3, 4] for our algorithms and 
w = [5, 10, 20, 30, 50] for Laplacian-based regularization. Such ranges were empirically found to 
be optimal regularization weight values for our data sets. Figure [3] shows the estimated transfor- 
mation RMSE for different values of w. Our adaptive regularization approach is accurate with an 
average transformation error below 1 voxel, whereas quadratic regularization has less satisfactory 
performance. 




(a) reference (b) source (c) result 



(c) difference before (d) difference after 



Figure 2: 3D non-rigid registration example. We register the source image (b) onto the reference 
image (a) to obtain the registered image (c). The difference image (d) between the reference and 
registered image is almost zero, which demonstrates the accuracy of the registration. The estimated 
average transformation error is ermse — 0.82. 
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Figure 3: Comparison of the non-rigid image registration performances of our algorithm and Lapla- 
cian (curvature) based method. We conducted 20 experiments for each of the regularization weight 
values. We used w — [0.5, 1.2, 2, 3, 4] for our algorithms and w = [5, 10, 20, 30, 50] for Laplacian- 
based regularization. The Laplacian-based regularization results were rescaled into the common 
range [0..5] for simpler visualization. Our adaptive regularization approach is accurate with aver- 
age transformation error below 1 voxel, whereas standard Laplacian-based regularization has less 
satisfactory performance. 



5 Related work 



Analogy with Kernel matrix completion: Tsuda et al. |fT9l proposed a method to complete a ker- 
nel matrix with auxiliary data. In their work, the input data is available only for a subset of samples. 



7 



and the kernel matrix derived from such data has missing entries. To complete the kernel matrix, 
the authors use an auxiliary kernel matrix derived from another information source. They minimize 
the KuUback-Leibler (KL) divergence between two kernel matrices and make use of the Riemann- 
ian information geometry, where the KL-divergence is defined by relating the kernel matrix to a 
CO variance matrix of Gaussian distribution. The KL-divergence allows to use the em algorithm [201 
(different from the EM algorithm of Dempster et al. [21]). In these terms, by minimizing KL{p\ \q) 
(or in matrix notations K L{ll^^\\Vl^^) in Eq.|9l), we are finding (on the manifold) a positive def- 
inite matrix that is closest to the given matrix Q,^^. Relating our formulation to the kernel 
completion one, we are simultaneously completing the inverse covariance matrix from the cur- 
rent (estimated) observations u* and the auxiliary matrix Q,^^. This allows flexibility on E^^, in 
contrast to its fixed form formulation. 

Analogy with ii-norm regularization: Consider a simple case of weighted L2 norm regulariza- 
tion. 

vamEUi, a) -L>(0, u) + V a,uj (24) 
w ^-^ 

i 

where > are unknown weights. This is equivalent to the assumption of a strict diagonal form of 
the covariance matrix S. To see the analogy, the last term here can be written as ^ Ojuf = u-^E^^u, 
where has elements ai along its diagonal. We assume the model distribution (7 to be a Gaussian 
with isotropic diagonal covariance of all ones, that is ^(u) oc e^H"" Following our derivation 
we can analytically solve for a^, ehminate it from the equation and achieve the following objective 
function 

min£;(u) = -D(Cl,u) + VluJ (25) 
w ^-^ 

i 

which is the Li regularized problem. Thus, optimizing the error function with L2 weighted norm 
(with unknown weights) is equivalent to the Li regularized problem (also called Lasso in regression 
problems |22|). This provides another interpretation of Li norm regularizers that recently attracted 
a lot of interests in the machine learning community ||23ll24ll . 



Analogy with Sparseness and Compression: In our final objective function (Eq. [T5]i, the last 
term represents the Li norm of the DCT (or FFT) coefficients of u. Li norm has 

been popular to enforce sparseness of the coefficients |25|. From this perspective, the standard 
quadratic regularizer can be seen as the one that penalizes L2 norm of the DCT coefficients, e.g. 
II Au||^ = u-^QK^Qu ~ ^ fc|(q^u)^ (compare to our ^ fcilqj'ul). Thus using our regularization 
approach we are also enforcing sparseness of the DCT coefficients of u. Sparseness of the estimated 
signal often leads to its better generalization properties [23J. As a few DCT coefficients include 
most of the signal information, sparseness of the DCT coefficient also forces higher compression of 
the estimated signal. 



Analogy with Adaptive Filtering: Finally, we draw the analogy of our approach with adaptive 
filtering. Optimization with a standard quadratic regularizer is equivalent to unregularized opti- 
mization followed by filtering, where the filter depends on the regularization operator 1261 . Indeed, 
Eq. I20I represents a gradient descent step, followed by the fixed low -pass filter (I + 7wA)^^ in 
frequency domain. In our adaptive regularization approach, the filter becomes signal-dependent: 

(d IQ'^ ' u'l+JmK) ^^^^ EU- Such a filter also resembles the Wiener filter, where the power spec- 
trum of the actual signal is |Q^u* | (signal from the previous iteration) and the noise power spectrum 
is 7i(;K. The important fact here is that the filter is adaptive, it changes at every iteration, whereas 
for the standard regularization operator, the filter is fixed. 



6 Discussion and Conclusion 

We introduced adaptive regularization approach. Instead of assuming a fixed regularization opera- 
tor (or a fixed prior distribution on the function/parameters), we estimate it. We assume the prior 
distribution on parameters to be close to the given model distribution in terms of KL-divergence. 
We constrain the prior distribution to be a Gauss-Markov random field, which allows us to solve 
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for the prior distribution analytically and eliminate it from the equation. The final objective func- 
tion appears to have a regularization term that penalizes the absolute value of u after applying an 
orthogonal transformation (Fourier or DCT). DCT approximates the optimal (in the decorrelation 
sense) Karhunen-Loeve transform with certain Markovian assumptions |8, 7|. Penalizing the abso- 
lute value of the decorelated vector u, we are also enforcing sparseness on u in terms of the basis 
functions (DCT, DFT), which leads to better compression and generalization properties. We also 
proposed the fast optimization algorithm with complexity of 0{N log N). 

Using our regularization approach, we achieved accurate non-rigid image registration results. Our 
method recovered challenging non-rigid transformation fields, whereas standard variational methods 
had less satisfactory results and usually converged to a "bad" local minimum. We still have to 
choose the regularization weight and the regularization operator (covariance matrix of the model 
distribution) similar to the standard quadratic regularization. Nevertheless, our approach does not 
constrain the regularization operator to a fixed form and allows it to be flexible. 
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